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Abstract. Wc develop a family of expanded mixed Multiscale Finite Element Methods (Ms- 
FEMs) and their hybridizations for second-order elliptic equations. This formulation expands the 
standard mixed Multiscale Finite Element formulation in the sense that four unknowns (hybrid 
formulation) are solved simultaneously: pressure, gradient of pressure, velocity, and Lagrange mul- 
tipliers. We use multiscale basis functions for both the velocity and the gradient of pressure. In the 
expanded mixed MsFEM framework, we consider both separable and non-separable spatial scales. 
Specifically, we analyze the methods in three categories: periodic separable scales, G— convergent 
separable scales, and a continuum of scales. When there is no scale separation, using some global 
information can significantly improve the accuracy of the expanded mixed MsFEMs. We present a 
rigorous convergence analysis of these methods that includes both conforming and nonconforming 
formulations. Numerical results are presented for various multiscale models of flow in porous media 
with shale barriers that illustrate the offlcacy of the proposed family of expanded mixed MsFEMs. 

Key words, expanded mixed multiscale flnite clement methods, non-separable scales, hybridiza- 
tion, two-phase flows 

AMS subject classifications. 65N99, 65N30, 34E13 

1. Introduction. There are many fundamental and practical problems involv- 
ing a wide range of length scales. Typical examples include highly heterogeneous 
porous media and composite materials with fine micro-structures. In practice, these 
scales are too fine to treat with direct numerical simulation as the computational cost 
far exceeds the capabilities of computers for the foreseeable future. As a result, it 
is a significant challenge, both theoretically and numerically, to treat problems with 
multiple-scales effectively. A variety of numerical algorithms, varying from upscal- 
ing to multiscale methods, have been developed to capture the influence of fine-scale 
spatial heterogeneity on coarse-scale properties of the solution. In most upscaling 
methods the coarse-scale model is developed by numerically homogenizing (or av- 
eraging) the parameters of the fine-scale model, such as permeability. The form of 
the coarse-scale model is typically assumed to be the same as the fine-scale model. 
Simulations are performed using the coarse-scale model, and coarse-scale quantities 
of interest are readily computed. However, with this parameter upscaling, fine-scale 
properties of the solution can no longer be recovered. In contrast, multiscale methods 
carry fine-scale information throughout the simulation, and the coarse-scale equations 
are generally not expressed analytically, but rather formed and solved numerically. 

Various numerical multiscale approaches have been developed during the past 
decade. A Multiscale Finite Element Method (MsFEM) was introduced in [23] and 
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takes its origin from the pioneering work j5] . Its main idea is to incorporate fine-scale 
information into the finite element basis functions and capture their effect on coarse 
scales via finite element formulations. In many cases, the multiscale basis functions 
can be pre-computed and used repeatedly in subsequent computations with different 
source terms, boundary conditions and even slightly different coefficients. In some 
situations the bases can be updated adaptively. This leads to a large computational 
saving in upscaling multi-phase flows where the pressure equation needs to be solved 
many times dynamically. There are a number of other multiscale numerical methods 
with different approaches, such as residual free bubbles |II] . the variational multiscale 
method [24] , two-scale conservative subgrid upscaling [4] and the heterogeneous mul- 
tiscale method [T^- Arbogast et al. [B] used a domain decomposition approach and 
variational mixed formulation to develop a multiscale mortar mixed MsFEM. Jenny 
et al. [15] have used a set of multiscale basis functions similar to [53] to develop a 
Multiscale Finite Volume Method (MsFV). A Multilevel Multiscale Method [IB] was 
proposed in the framework of the Mimetic Finite Difference Method that efficiently 
evolves a hierarchy of coarse-scale models. In recent years, multiscale methods (e.g., 
[ini [m dH] ) have been developed to address high-contrast in multiscale coefficients. 

The mixed multiscale method was first developed by Arbogast et.al. [3] as a 
locally conservative subgrid upscaling method, and was later analyzed and extended 
in [S] using a variational multiscale formulation. Chen and Hou developed a local 
multiscale basis equation for velocity and combined it with a mixed finite element 
formulation to propose a mixed MsFEM [M]. The mixed MsFEMs retain local con- 
servation of mass and have been found to be useful for solving flow equations in 
heterogeneous porous media and other applications. However, in many practical sit- 
uations, the permeability may be very low, and even vanish in some regions of the 
domain (e.g., permeability in shale). In this case its inverse is not readily available for 
use in the standard mixed MsFEM formulation, and hence, the direct application of 
these methods to practical problems may be restricted. To overcome this limitation, 
we propose a family of expanded mixed MsFEMs for these multiscale applications. 

Expanded Mixed Finite Element Methods (MFEMs) were studied in classic finite 
element spaces in the past. Chen [T3] developed and analyzed expanded MFEMs for 
second-elliptic equations, and obtained optimal error estimates for linear elliptic equa- 
tions and certain nonlinear equations in standard finite element spaces. Woodward 
et al. [32] performed an error analysis of an expanded MFEM using the lowest-order 
Raviart-Tliomas space for Richards equation. In [^ Arbogast et al. established the 
connection between the expanded MFEM and a certain cell-centered finite difference 
methods. A dual-dual formulation was introduced in [3T] for the expanded MFEM 
with Raviart-Thomas spaces. The expanded formulation extends the standard mixed 
formulation in the sense that three variables are explicitly approximated, namely, the 
scalar unknown (e.g., pressure), its gradient and the velocity. The expanded MFEMs 
are suitable for the cases where the tensor coefficient of the underlying partial dif- 
ferential equations is small and even partially vanishing inside cells, which may be 
viewed as an extreme case of high-contrast of coefficients. 

Porous media formations are often created by complex geological processes and 
may contain materials with a widely varying ability to transmit fluids. The permeabil- 
ity of the porous media may change dramatically, from almost impermeable barriers to 
highly permeable channels. This high-contrast of the permeability creates significant 
challenges for the simulation of fiow through porous media. If the simulation fails to 
capture the influence of these barriers (e.g., shales) and channels (e.g., sand lenses). 
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errors in the quantities of interest will likely be unacceptably large. Moreover, if the 
low permeability regions are thin with a small inter-barrier spacing, a fully resolved 
fine-grid discretization may be difficult to obtain and costly to apply. Thus, a multi- 
scale simulation on a coarse grid may be necessary for these situations. To simulate 
the flows on a coarse grid when the permeability is almost vanishing in some regions, 
a reduced-contrast approximation was developed in MsFEM [16] to lower the variance 
of the coefficients. In [2], Aarnes et al. proposed an automatic approach to detect the 
barrier and iteratively split the coarse cells around barriers to obtain improved solu- 
tions in mixed MsFEM. For these barrier situations, the proposed expanded mixed 
MsFEM can provide an easy and direct approach to perform the computation. 

The purpose of this work is to develop a framework of expanded mixed MsFEMs 
and rigorously analyze the convergence for different multiscale cases. The work en- 
riches the studies on mixed MsFEMs from previous works [TJ [21 [M] . The multiscale 
phenomena can be roughly classified into two categories: separable and non-separable 
spatial scales. In the case of separable scales one can localize the computation of mul- 
tiscale basis functions. However, these approaches usually produce resonance errors 
Strategies exist to reduce the resonance errors, such as the oversampling technique 
introduced in |23] . Recently, a new technique was proposed in [22) . where a zero-th 
order term is artificially added to the standard multiscale basis equation [53] to make 
the associated Green's function decay exponentially, and consequently, the resonance 
error can be reduced. However, these strategies for reducing the resonance errors 
usually result in a nonconforming FEM, and moreover, they do not remedy the poor 
performance observed when the there is no scale separation. Instead, when media 
exhibit strong non-separable scales, some global information is needed for represent- 
ing non-local effects. If the global information captures all relevant information from 
different scales, then resonance errors are removed and approximation accuracy is sig- 
nificantly improved [TJ [TOj [30] . If global information (or field) is used, we refer to the 
multiscale methods as global multiscale methods; otherwise, we refer to them as local 
multiscale methods. We consider the expanded mixed MsFEM for both separable 
and non-separable scales. Multiscale basis functions are employed for both the veloc- 
ity and the gradient of the scalar. We present convergence analysis for three typical 
multiscale cases: periodic highly-oscillatory separable scales, G— convergent separa- 
ble scales, and a continuum of scales. For the expanded mixed MsFEM, we consider 
both conforming multiscale bases and nonconforming/oversampling multiscale bases. 
The computation of expanded mixed MsFEM is very similar to the standard mixed 
MsFEM, but the expanded mixed MsFEM gives more information about the solution 
and provides an accurate approximation of the scalar field's gradient as well. This 
gradient is often needed in applications, for example, the pressure gradient is used for 
flows with gravity in computing upstream directions. 

We consider the proposed expanded MsFEMs in the non-hybrid (standard) form 
and the corresponding hybrid form. The hybridization of the mixed formulation is 
of increasing importance to the study of mixed methods. Hybridization was initially 
devised by Fraejis de Veubeke jl9| as an efficient implementation technique for mixed 
finite elements. Specifically, hybridization localizes the mass matrix on each cell, and 
hence, enables the local elimination of the velocity to obtain a sparse positive defi- 
nite system. Later hybridization was applied to produce a better approximation of 
the scalar unknown over each cell. For example, Arnold and Brezzi [8] showed that 
using the Lagrange multiplier unknowns introduced by hybridization in a local post- 
processing procedure improved the accuracy of the scalar unknown. In the framework 



4 



L. Jiang, D. Copoland and J. D. Moulton 



of expanded mixed MsFEMs, we first show the equivalence between each expanded 
mixed MsFEM and its hybridization, and then analyze the convergence of these meth- 
ods. In particular, we obtain the convergence rates of the Lagrange multipliers for 
each of the three multiscale cases. 

The rest of the paper is organized as follows. Section 2 is devoted to formu- 
lating the expanded mixed formulation for a model elliptic equation in an abstract 
framework. In Section 3, we present multiscale basis functions and formulations of 
expanded mixed MsFEM for separable scales and continuum scales. In Section 4, we 
present convergence analysis for both expanded mixed MsFEMs and their hybridiza- 
tions. Local multiscale methods and global multiscale methods are analyzed for the 
expanded mixed MsFEMs. In Section 5, the expanded mixed MsFEMs are applied 
to various multiscale models of flow in porous media to demonstrate their efficacy. 
Finally, some comments and conclusions are made. 

2. Background. In this section we highlight the formulation of the expanded 
MFEM, and summarize important notation. We focus on the dual-dual formulation 
of the continuum problem because it most naturally supports our analysis. We review 
important results of the corresponding discretization as well. 

2.1. The continuum expanded mixed formulation. Let 17 be a bounded 
domain in M.'^, d ^ 2,3, with Lipschitz boundary dft. For a subdomain D, m > 
and 1 < p < oo, W™'P{D) and LP{D) denote the usual Sobolev space and Legesgue 
space, respectively. The norm and seminorm of W™'^ are denoted by || • ||m.p.n and 
I • Im.p.D, respectively. When p = 2, W"^'^{D) is written as H"^{D) with norm || • \\m.D 
and seminorm | • \m.D- The norm of L'^{D) is denoted by || • ||o.d- We also use the 
following spaces 



For a normed space X on f2, || • \\x denotes the underlying norm on X. We denote 
X{D) to be the restriction of X to subdomain D. In the paper, (•,•) is the usual 
L'^iP) inner product, and (•, ■)d is the L'^{D) inner product. 
We consider the following elliptic problem 



We may assume k is symmetric positive definite. Let 9 = Vp and u = — fc0. The 
equation (j2.ip can be rewritten as 



Here u usually refers to the fluid velocity. We define the following function spaces for 
solutions: 



Hh{D) 
H{div,D) 
H°{div,D) 



{q e H\D) : q\aD ^ 0} 

{v e [L^{D)f : \\v\\mv^d ■■= \\v\\o.D + \\div{v)\\o^D < oo} 
{v e H{div, D) : V ■ n\QD = 0}. 




(2.1) 




(2.2) 



Xi = [L^in)f, X2 = H{div, n), X ^Xix X2, Q = L^{n). 
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The expanded mixed formulation of (|2.ip reads: find {9,u,p) in Xi x X2 x Q, satisfying 

(2.3) 



{9, t) + (p, divT) = Vr e X2 
[divu, q) = (/, g) q^Q 



We have the following theorem, which addresses the relation between the solution of 
(1231) and the solution of 

Theorem 2.1. (See Theorem 3.5 in If {{B,u),p) € X x Q is the so- 

lution of l\2.S\) . then p G i?o(f^) is the solution of i2.1\) . satisfying 9 = Vp and 
u = ~k9. Conversely, if p £ Hq{VI) is the solution of i2.1\) . then h2.S\) has the 
solution (^{0,u),p^ € X X Q satisfying 9 — Vp and u = —k9. 

For theoretical analysis, we introduce the operators associated with (|2.3p . Let 
Ai : Xi — > X[ and Bi : Xi ^ 
space. These arc defined, respectively, by 



where the notation ' denotes the dual of a 



[A,{9),£]:={k9,0, VeeXi, 
[B^{9),v]:=(9,v), yveX2, 

from which we define A : X — > X' by 

[A{9,u),{i,v)] [A,{9),i] + \Bl{u),^] + [B,{9),v], 

where B^ is the transpose operator of Bi. Let O denote the null operator. Then A 
can be rewritten as 



A = 



Ai Bl 
Bi O 



Let B : X2 — > Q' and F £ Q' . They are defined, respectively, by 

[B{u),q\ {divu,q), \/q E Q 
[F,q] (/,q),VgeQ. 

Using this operator notation, the expanded mixed formulation (|2.3p is equivalent to 
finding (^{9,u),p^ in X x Q such that 

O] r 61 ] I O' 

(2.4) 

Equation (|2.4p is called a dual-dual mixed formulation of (|2.3p PO] , since the operator 
A itself has a similar dual-type structure. Because the operator B satisfies a continu- 
ous inf-sup condition |12j and the operator Bi satisfies a continuous inf-sup condition 
on kcr(i3), the system (|2.4p has a unique solution (^{9,u),p^ in X x Q. Moreover, 
there exists a C > 0, independent of the solution, such that [21] 



Ai 


Bl 


0' 




' 9 ' 




' ' 


Bl 





B* 




u 










B 







P 




F 



\{{0,u),p)\\xy<Q<C\\f\\o,n.. 



(2.5) 



where 



\\{{9,u),p)\\xxQ II^'ILyi + \\u\\x2 + WpWq ■= M\o,n + \\u\\div.n + \\p\\o,i 
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2.2. The discrete expanded mixed formulation. Let Aij,, Bi ^ and be 
approximations of operators Ai, Bi and B, respectively. Let F}i be an approximation 
of F. Then the numerical formulation of (|2.4p reads 



Aim 
Bijt 
O 



B 



i./i 
O 

Bh 



O 
O 





' Oh ' 

























(2.6) 



Let [•, •]/! be a numerical inner product. Let Xi^, X2.h and Qt be finite dimen- 
sional approximation of Xi, X2 and Q^, respectively. Set Xh Xi_h x X2,h- Then 
we have the following abstract result for well-posedness of (|2.6p . 

Lemma 2.2. Assume that 
(1) there exists a positive constant Ci independent of h such that for any G Qh, 

[Bhivh),qh]h 



sup 

Vh&X2.h\{0} 



>Ci||<7h|lQ; 



(2.7) 



\'"h\\X2 

(2) there exists a positive constant C2 independent of h such that for any Vh G kcr(_B;j) 

[Bi.h{S.h),Vh]h 



sup 



\ih\ 



Xi 



(2.8) 



(3) there exist positive constant C3 and constant C4 independent of h such that 
CsWOhW'x, < [AiM{0h),0h]h < CiphWj,^ for G X^.,,. 



Then there exists a unique {OhTUf^Ph) £ Xi^h x X2,h x Qh solution of 
We note that || • \\x2 should be an element broken norm, provided that X2.h is 

a nonconforming space in X2, i.e, X2,h 't- ^2- We have the following Strang- type 

lemma for the problem (|2.6p . 

Lemma 2.3. flQj / Let ({9,u),p) e X x Q and ({9h,Uh),ph) G Xh x Qh be the 

unique solutions of {2.4^ and h2. 0fl . respectively. Then there exists C = C(||i?i||, Ci, C2, C3, C4) > 

0, where Ci, C2, C3 and C4 are defined in Lemma \2.2l such that 



(2.9) 



{{d,u),p) - {{Oh,Uh),Ph)\\xxQ < 

^ \\{iO,u),p) 
Bi^hiO) - Bl{p),vh]h 



C{ inf \\{{0,u),p)~{{^h,vh),qh)\\xxc 

{{ih.Vh).qh)eXhXQh 



+ sup 

Vh&X2,h\{0} 



\\Vh\\x2 



}• 



Again, || • jjxa in Lemma 12.31 should be an element broken norm, provided that 
X2,h is a nonconforming space in X2. The second term on the right hand side of (|2.9p 
is the so-called consistency error. It is easy to check that 

[-Bu,{e)-Bl{p),Vh]h _^ _ [-B,{0),ih] 



sup 

VheX2,h\{o} 



sup 

5hekcr(B,0\{0} 



\Vh\\X2 

Consequently, if kcr{Bh) C ker_B, then 

[-B,,h{0)-Blip),vh]h 



Vh X2 



sup 

Vh&X2^h\{0} 



\vh\\x2 



0, 



which is the case for a conforming expanded mixed FEM. 

As a general remark, the generic constant C is assumed to be independent of the 
mesh diameter h throughout the paper. 
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3. Formulations of expanded mixed MsFEMs. In this section, wc describe 
how to construct multiscale basis functions for expanded mixed MsFEMs. 

Let %h be a quasi-uniform partition of £7 and if be a representative coarse mesh 
with diam(K) = hx- Let h = max{hK, K € T/i}. Let (Vh,Qh) be a classic mixed fi- 
nite element space pair such as the Raviart-Thomas, Brezzi-Douglas-Marini or Brezzi- 
Douglas-Fortin-Marini spaces (cf. [E]). In this section, wc present multiscale basis 
functions for standard local expanded mixed MsFEM, oversampling expanded mixed 
MsFEM, and global expanded mixed MsFEM. 

3.1. Local expanded mixed MsFEM. Let Vh{K) := VhIk (the finite element 
velocity space localized on K) and G Vh {K) be a velocity basis function defined 
on K. Following the mixed MsFEMs developed in [T3], we define the multiscale basis 
equation for the local expanded mixed MsFEM in the following way: 



+ kr)^ =0 in K 

* ' — in r\ 

(3.1) 



r;^ - V(^^ = in K 
div{'>p^) = divx^ in K 
■ riQK = x'^ ■ ndK on dK, 

where ngK is the outward unit normal to dK. We can obtain different mixed MsFEM 
basis functions corresponding to different choices of x^ ■ Chen and Hou [13] choose 
X^ to be the lowest Raviart-Thomas basis. The weak expanded mixed formulation 
of dnil) is to find ((77^, V'^), 0^) e X{K) X Q{K)/R such that 

(?7^, r) + {(j>^, divT) =0, Vr e H°{div, K) (3.2) 
{div^l;^,q) = {divx'',q), yqeL^K) 

with ijj^ ■ UQK = X^ ■ ^dK on dK. The function 77^ is a basis function for the gradient 
variable 9 and -0^ is a basis function for velocity u. We define finite element spaces 
for the local expanded mixed MsFEM as follows: 

Xlh — {0 e [L^in)f : 0\k e span{r]i^} for each K e T^} 
X2,h {" •= H{div, fl) : u\k G span{^^} for each K e 1h} 
Xh X[ ,^ X X2 ,^. 

The finite element space for pressure in the expanded mixed MsFEM is Qh, which 
is the same as in the classic mixed finite element pairs. On each coarse element 
K, the velocity basis functions are associated with its faces. Therefore, velocity basis 
functions of X2 have the support of two coarse elements sharing a common interface. 
However, the basis functions of X\ ^ are supported in only one coarse element. Each 
basis function of X\ ^ is associated with a velocity basis function via ip^ = —kr]^, 
but it is restricted to only one coarse element K. 

The local expanded mixed MsFEM formulation of the global problem (|2.3p reads: 
find ((6'^,M'h),p'h) e Xl X Qh such that 

(fceLa) + «,a) = V4 

{elrh) + {pldivrh)^0 yrn^Xlh (3-3) 
{divuh, qji) = (/, Qh) Q ^ Qh 
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subject to the global boundary conditions. 

It is well known [13) that the linear system of equations in the expanded mixed 
FEM produces an indefinite matrix, which is a considerable source of trouble when 
solving the linear system. We can introduce Lagrange multipliers on interfaces of 
cells and localize mass matrix to each element to obtain a sparse symmetric positive 
definite system by elimination, which is suitable for many linear solvers. This will 
give rise to a hybrid formulation. To this end we require some further notations. 

For the hybrid formulation we define the MsFEM velocity space to be ^, 
whose normal components are not necessarily continuous on J^^, the collection interior 
interfaces of 'Xh- They are defined as follows: 



{u e [L'^{fl)]'^ : u\k e span{V'5} for each K £ X/i} 



We note that basis functions in ^ do not require normal continuity. The basis 
function in X2 is supported on a coarse grid, and this is different from the basis 
function in X2 whose support is two coarse cells sharing the common face. We note 
that 



X 



2,h 



xl„nH{dtv,n). 



The finite element space for the Lagrange multiplier is defined as 



Define the operator C'^ : Xj ,j — > 



for each e £ Fh, 7r|e — for e C 9ri}. 



(n',)' by 



[ci{TH)M 



K 



{Th ■ nOK, T^h)dK, e X2 jj, TTh e II/j 



The hybridization of local expanded mixed MsFEM formulation of 



p'^, A;J e Xf^ X X^ f^ X Qh X n',^ such that 



K 



{pi,divTn)^[Ci{Th),A] ^Ti^eXl^ 

{f,qh) q^Qh 



qh) 



K 



(3.4) 



reads: find 



(3.5) 



3.2. Oversampling expanded mixed MsFEM. If the equation (j3.ip is solved 
in a block S larger than K , and the interior information of the solution is taken to 
construct the basis functions in K, then this results in the oversampling technique 
introduced in [23l[T4). We note that oversampling MsFEM is a modified local MsFEM. 
Using the oversampling can reduce resonance error. 

Here we follow the outline in |14| to present the oversampling multiscale basis 
functions. Let ((77^, V'^), (/)^) e X{S) x Q{S)/R be the solution of the equation 



4 



k4 = 







m 
in 



S 
S 

div{ij}^) — divx^ in S 
■ ngx = ■ ridK on dS. 



(3.6) 
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Let j and / be indices of faces of K. We define 

where the constants Cji are chosen such that 

xf = Y.'^^l^l\^■ 



(3.8) 



We define 



^2% ■■= {u e [L%n)f : u\k e span{V'^'} for each K e 



We introduce an operator Ai h ■ — ^ whose local form is defined by 

MH\K{Y.c,i;^J^Y.^,xf. (3.9) 



Then we define finite element spaces for the oversampling expanded mixed MsFEM 
as follows: 

X?% ■■= {0 e [L^{n)f : 0\k e span{fj^} for each K e T,,} 
^2% {uh e : MhUh £ 14 C H{div, n)}. 

Here, AitUh G Vh is to impose some intrinsic continuity of the normal components 
of velocity multiscale basis functions across the interfaces Fh- However, X!^^ <f_ 
H{div, n) and so the oversampling multiscale method is nonconforming. The ex- 
panded oversampling mixed MsFEM formulation of the global problem (|2.3p reads: 
find (6'^^^t^^p^/) e X°^^ X X°-^^ x Qh such that 



K 

^{divu°f,qh)K = (/, Qh) Qh e Qh 



(3.10) 



subject to the global boundary conditions. 

We define operator C^" : X°% — > (H'^)' by 

[Cr(T/.),7rJ := ^(A^/.(r^) • n,7rh)sK, Vr^ G la"!,^'. € H'^. 

K 

Then the hybridization of ((3J0| is to find {ei^,ul',pf, XI') E XfJ, x X"'), x Q,, x 
such that 

(fcC,60 + (^ir,a) = o vae^a 



K 

'^[divu°f^' ,qh)K = (./, gft) qh e 



(3.11) 



if 
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3.3. Global expanded mixed MsFEM. To remove the resonance error and 
substantially improve the accuracy, we can use global information to construct the 
multiscale basis functions. Suppose that there exist global fields ui,--- , u^r that 
capture the non-local features of the solution of equation (|2.ip . Following the global 
mixed MsFEM framework proposed in [I], we define the multiscale basis equations 
for the global expanded mixed MsFEM by 



in K 
in K 

K 



div{'ipi = divx in K 



on dK, 



(3.12) 



where bf^ 



IdK -naxds 



Ui ■ UQK and i = 1, N are indices of the global infor- 



laK ' "^SKds 

mation. References [J [101130] provide some options for the global fields. To reduce the 
computational cost of using global fields, we can pre-compute these fields ui, • • • , 
at some intermediate scale [26] . 

The finite element spaces for the global expanded mixed MsFEM are defined by 

^f.h '■= e [L^{n)f : 0\k e span{r]^^} for each K € S,,} 



Y9 



{u e H{div, n) : G spanlipf } for each K e 1^} 



Consequently the global expanded mixed MsFEM formulation of the global problem 
(lO) reads: find ((6*^,0,^^) € x Qh such that 

(fc^^^,a) + K,a) = o vaexf, 

{9lTh) + {pldivTh) = yrh^Xl^ (3.13) 

{divul,qh) ^ {f,qh) q G Qh 

subject to the global boundary conditions. 

We define the MsFEM velocity space h ^-"^"^ ^^^"^ finite element space of the 
Lagrange multipliers for the hybrid formulation of p.l3p by 

X^ , -.^ {u e [L^in)f : u\k e spaniijf;^} for each K G T,,} 



Ul := {t: e L''{J^h) : 7:\, G X. 



2.h 



for each e G J^i, 'I'le = for e C dil}. 



The functions in X| may not have normal continuity. The relation between X| ^ 
and X| ^ is expressed by the identity X| ^ = X| ,^ n H{div, fi). 



Let the operator : X. 



2,h 



(n^)' be defined in a similar way to (|3.4p . Then 
the hybridization of the global expanded mixed MsFEM formulation of p.l3p reads: 
find {91, ul,plXl) e X^^x Xl,^ X Qh x Hf^ such that 



(fc0l,e^) + (u^,a) = O vaex,% 



a: 

E 



{divul, qh) = (/, g/0 q £ Q/i 



(3.14) 



[C,l(i2^),7r„]=0 n^eUl 
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We note that 



(|3.11l) and p.l4[) can be rewritten as the operator formulation 















' h ' 




o 


Bi^h 





Bl 


-CI 








o 


O 


Bh 










Ph 




Fh 


O 


-Ch 










. _ 




o 



(3.15) 



This is actually the hybrid approximation of the operator equation (|2.6p . 

4. Convergence analysis. We will focus on the analysis of the local expanded 
mixed MsFEM. The analysis of the oversampling expanded mixed MsFEM and the 
global expanded mixed MsFEM is very similar to the local expanded mixed MsFEM 
but requires additional notation and definitions. We sketch the analysis of the over- 
sampling expanded mixed MsFEM and global expanded mixed MsFEM and present 
their convergence results. For the local expanded mixed MsFEM, we consider micro- 
scales that can be treated in terms of periodic homogenization and G— convergent 
homogenization. 

For the analysis, we will focus on the case when the velocity space on Vh{K)^ 
is the lowest-order Raviart-Thomas space (i?7o) and ^ Vh{K). We use to 
build the source term and boundary conditions for the MsFEM basis defined in equa- 
tion (O). 



4.1. Inf-sup condition for expanded mixed MsFEMs. In this subsection, 
we discuss the inf-sup condition associated with problem p.3p . To simplify the pre- 
sentation, we define three discrete operators associated with p.3|) . Specifically, let 



are defined as 



X[j^, Bi^h ■ Xj^j^ — 

[BiM{e),v] 
[Bh{u),q] 



X' and Bh : X^,^ 



Q'l^, then the operators 



--{e,v), V«eX^,,, 

-- {div{u),q), \fqeQh- 

By Lemma 12.21 we require that the operator Ai^h is bounded and positive for well- 
posedness of problem p.3p . If k satisfies the assumption 



fcminlCI' < C^k{x)C < fcmaxICI' for VC G K'^ and V.T e n, 



(4.1) 



then the operator Ai^h is bounded and positive. However, we can weaken the assump- 
tion (|4.ip : for example, k{x) may be partially vanishing inside fine cells and positive 
and bounded everywhere else. Then the operator Ai h is still positive and bounded. 
For simplicity of presentation, we use assumption (|4.ip for the analysis. 

We have the following lemma for the discrete inf-sup condition of problem p.3p . 

Lemma 4.1. For any qh € Qh, there exists a positive constant Ci independent of 
h such that 



[Bh{vh),qh] - „ II I 

sup r r > Ci\\qh\ 

VHexl^\{a} \\vh\\x2 



(4.2) 



Let kcT(Bh) be the kernel of the operator Bh- Then for any Vh G ker(i?/i), there exists 
a positive constant C2 independent of h such that 



sup 

€heA-}^\{0} 



[Bi.h{^h).Vh] 



\ih\ 



>C: 



Xi 



2\\Vh\\x,- 



(4.3) 
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Proof. Wc define € V/i, whose localization on K is v^^Ik — o-fxf € Vh[K), 
where j is the index of a face of K. Consider a map M whose local form is defined 
by MkvI = J2j o.f'^xo ^ ^2,hi^)- Then M : Vh — > Xi,^^ is a one-to-one map. It is 
easy to check for any v*^ £ Vh, 

div{Mvl) = div{vl) in K. (4.4) 

Let z-^' = «f '/'x, • Then Mkv^ = -kV and we have 

WkK\\Ik = I k-^z^ -k^/z^dx <C I kS/z^ ■\/z'^dx 

J K Jk 

= ~C ( MkvI ■ Vz^dx = C( /" div{MKvl)z'^ dx - ( [MkvD ■ nz'^ds) 

Jk Jk Jok (45) 

= C( / div{vl)z^dx - vl- nz^ds) =~C < • Vz^dx 
Jk JdK Jk 



C vl- k-HlKvldx < C\\vt\\o,KmKvt\\o. 



k- 



k 



This gives |jMifU^||o,if < C|k;*Jlo-A'- Consequently, combining with (|4.4I) implies that 
for any < e V^*, 

\\Mvl\\x,<C\\vl\\x,. (4.6) 

Consequently, since the inf-sup condition holds for the classic mixed FE pair (V^, Qh), 
it follows that for any qh € Qh 

{divvh,qh) ^ {div{Mvl),qh) 
sup — — > sup — II ^ ^ ' 

v^eX'MO} \\Vh\\X2 v'^eV^MO} P^^'hllXa 

(4.7) 

. 1 {divvl^qh) C* 

^v*£Vh\{0} \\Vh\\X2 ^ 



where we have used (|4.6p in the second step. The proof of (|4.2p is complete. 
It is obvious that by definition 

ker(Sft) = {vh e X^j, : div{vh) 0}, 

and hence kei{Bh) C C X[ ^. Then for any Vh G ker(_B;i), 

[Bi.h{^h),Vh] {£.h,Vh) II II II II 

™P nT~n = 7m — = \\vh\\o.n = \\vh\\x2, 



which completes the proof of (|4.3p . □ 

Using Lemma 12.21 Lemma 14.11 and Lemma 12.31 gives the following theorem. 

Theorem 4.2. The problem \S.3\) has a unique solution ((^h., & Xl^y. 
Qh- Let (^{6,u),p^ be the solution of i2.3\) . Then there exists a positive constant C 
independent of h such that 

\\{{e,u),p)~{{ei,ui),pi)\\x.Q <C inf \\{{e,u),p)~{{ih,Vh),qh))\\x.Q- 

(4.8) 
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Following the proof of Lemma [4.11 and using the argument of Lemma 4.4 in 
we can obtain the inf-sup condition for the ovcrsampling expanded mixed MsFEM 
problem (|3.10p . Then applying Lemma [2731 to (|3.10p . we have the estimate 

Y,\mn),p) - ((c,o,pr)iUwxQw 

K 

-''^i,......)eW%^x,:,^Q.^^^^^ (4.9) 

_^ -{d,Vh) - Y.Adiv{vh),p) 



The last term of (|4.9p is the consistency error caused by the oversampling velocity 
space. 

Under suitable conditions described in [1], we can similarly obtain the inf-sup 
conditions for the global expanded mixed multiscale finite element problem (|3.13p 
and the Cca-type estimate, 

\\{{e,uipy{{eiui)y^)\\x.Q <c inf \\{{e,u),p)-{{in,vh),qh))\\x.Q. 

[{(.h,Vh),qh))eXly.Qh 

(4.10) 

4.2. Equivalence between expanded mixed MsFEM and its hybridiza- 
tion. In this subsection, we show that expanded mixed MsFEM formulation and its 
hybridization produce the same vector functions (i.e., the gradient of pressure and 
velocity). The result will be helpful to proceed with the convergence analysis for 
expanded mixed MsFEM and its hybrid formulation. For simplicity, we focus on the 
case of local mixed MsFEM. 

We first have the following lemma for operator C^. 

Lemma 4.3. Let C\ be defined in \3.4\l - Then 

kerC^ = Xi,„ ker(C7;j* = {0}. 



Proof. By definition of X\ ^ and (|3.4p . it immediately follows that 
Let K* be any cell in T/j and e* be any face of K* . Let r* G /j be such that 



and defined on K* by 



for e* 



where S VhiK*) is any basis function associated with face e*. Then [{C\^)*{'Kh), r^] 
means that (x^ • n, nhje' = for any basis function x^, G Vh{K*) associated with 
face e*. The classic mixed FEM theory implies that -Kh = on e*. Since K* is any 
cell in T/i and e* is an arbitrary interface, tt/i = on J-^. This completes the proof. 
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□ 

Next we show that problem l\'6.'6\i and problem (|3.5p produce the same solution. 
Theorem 4.4. Problem p.3p has a unique solution ((^'^, wjj,^'^) e x Qh- 
Problem \3.5]) has a unique solution (^^^ , , p'^ , A^) g X{ ^xXg ^xQ/iXll'^. Moreover, 

m,ui),pi) = {{elui)J^). 

Proof. 

Following the proof of Theorem 3.1 in [7] and using a standard duality argument, 
we obtain that 

||eLllo,n + ||tiLl|o,o + ||pLllo,n <Ci|/||o,o, 

which shows p.3p has a unique solution by setting / = 0. 

Let be the unique solution of (|3.3p . We define a linear operator by 

L(r,) = r,) + ^(p5„ MvTn)K, Vt„ € 1^^,. (4.11) 

K 

It is clear that L{Th) ~ for all Th E X2 ^, i.e, Xi, ^ C kerL and 

L e [ker(C[)]^ := £ (l^J' :< >= 0, Vr„ G ker(C^)}. 
By the Closed Range Theorem, it follows that 

[ker(C^r =5R((Ci)*). 
This means that there exists a A/i G II'^ such that {Cj^)*{Xh) = L, which gives 

i(r/.)-[C,Ur,),A$J Vr^el^,. (4.12) 
Then by Lemma [4.3) A'^ is unique. Combining (j4.1ip and ()4.12p gives 

(Si rn) + Y,{pldivTH)K - [Ci(r„), AJJ Vr„ G X^,,. (4.13) 
if 

The equation (|4.13p is exactly the same as the second equation in p.Sp . The fourth 
equation in (|3.5p and Lemma 14.31 show that Uh € ^2 h ■ uniqueness of solutions of 
p.3p and uniqueness of A'^, we immediately obtain that ((^^, ■s'^),p'^) is unique and 

□ 

By a procedure similar to Theorem 14.41 we have the following theorems for over- 
sampling expanded mixed MsFEM and global expanded mixed MsFEM, respectively. 

Theorem 4.5. Problem \3.11]) has a unique solution {91'^ ^ufi; ,p'^ .Xfi') £ X°'^f^ x 

X°% xQhX Moreover, ((07/,^^),^^) = {K' ,K').vl') , ^here {{91^ ,ul^),pl^) 
is the unique solution of problem iS.lOp . 

Theorem 4.6. Problem ^3.14^ has a unique solution (^^, A^) e X^ x 

Xl^ X X nl Moreover, {{9lul),pl) = {{elul),pl) , where {{9lui),pl) is the 
unique solution of (|3.13p . 

Thcorcm l4.41 Theorem l4.5l and Theorem l4.6l allow us to identify ((^^, "'JiP'^) with 
(K-O-fU' with ((0-,O,p-), and {i9lul),pl) with ((ef„<),p^), 
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respectively. So we can drop the upper bars in corresponding hybrid systems (|3.5p . 
p. lip , and (|3.14p . respectively. We note that uj^ ~ vJ-^^ is an identity in the sense of 
vector-valued functions. However, the corresponding coefficient arrays on their nu- 
merical representations are not equal to each other, because dimlX^ ^) 7^ dim^Xl^ y^). 
The same relationships hold for the oversampling expanded mixed MsFEM and the 
global expanded mixed MsFEM as well. Thus, to simplify the notation, we will drop 
the upper bars from p.5|) . p. lip and p.l4|) in the rest of the paper. 

4.3. Error analysis for the periodic case. In this subsection, we consider the 
case fc(x) = fc(f), where the parameter e characterizes small scales representing small 
physical lengths. To highlight the dependence on the small parameter e, we rewrite 
the solution of (|2.3p as ((6'c , Uf), p^) . We present an error analysis when the e scales are 
periodic. The local expanded mixed MsFEM and the oversampling expanded mixed 
MsFEM will be considered in the following subsections. 

We briefly recall the relevant homogenization results for the periodic case. Specif- 
ically, let y = ^ and k{y) be a periodic function in the unit cube Y = [0, 1]''. In this 
case, we can compute k* in the following way. Let N = {Ni,--- ^Nd} solve the 
auxiliary equations 

{-diVy{k{y)VNi) ^ diVy{k{y)ei) in Y 

\my))Y =0. ^ ■ ' 

Here (i = 1, • • • , d) is the unit vector in W^. Then the homogenized tensor k* is 
defined as (see, e.g., p7] ) 

k* ^ (fc(VAA + /))y. 

Let p* solve the homogenized equation of (|2.ip by 

in O 

(4.15) 




We define the homogenized velocity u* by u* —k*Vp*. 

Since we consider G Vh{K) := RTf){K), we have divx^ = ^ and x^'"e — ]^ 
on e and otherwise. For this case, Qu is a piecewise constant space. 

4.3.1. Analysis of local expanded MsFEM. Leti;^ : [H'^{Vl)Y — > Vh be the 
interpolation operator of RTq for velocity. Let Vh ■ i^(0) — > Qh be an orthogonal 

projection. We define the multiscale interpolation operator I/i ; [iJ^(O)]'* — > X. 
for velocity by 



I 

2,h 



where "0^^ is a multiscale velocity basis function associated with interface e. We define 
the interface projection Vd '■ L'^i^h) — > n'^ such that 

{Tap -p,nh)e^O Vtt/i e nj, and Ve G J'h- 

This interface projection implies that Vdp\e ~ iJeP'^^)Xe ■ n-e- We define a discrete 
norm on L-{J-h) by 



2 

KeleCdK 
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Lemma 4.7. If is a lowest Raviart- Thomas element, then the corresponding 
multiscale basis function tp^^ satisfies 



(4.16) 



Proof. By basis equation p.ip . it follows that 



1 



div{kV^^ ) = TT^ in K 



K\ 



-kVcl)^^ ■ Tie 



on e 



-k\7(j)^ -n^O ondK\e 



= 0. 



Multiplying (p-^ and using integration by parts for the above equation, we have 



1 



^1 



0. 



The Cauchy-Sehwarz inequality implies that 



(fcV0«,V0«) = - 



0,e- 



Applying the traee theorem and sealing argument |T2], it follows that 



|er^||^^J|o,e < C\e\-^hl\\y4>l\WK < C\e\-^hlUl\\o,K, 
Therefore, by (|4.17p and (|4.18p we have 



i.K I 



\,K- 



This eompletes the proof. □ 

We define the homogenization equation of the basis equation p.ip by 



r^r-^ - V0*'^ = 



= divx^ in K 



*,K 



noK = x'^ ■ n-dK on dK . 



By Theorem 12.11 and the uniqueness of solution of ()4.19p , we have 

Let (x) be the solution (up to a constant) of the equation 

i-div{k,{x)\'w^) ^div{Ilu*\K) in K 

\ -K{x)Vwf ■ nsK = {IIu*\k) ■ nan on dK. 



(4.17) 



(4.18) 



(4.19) 



(4.20) 



(4.21) 
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Then the homogenized equation of (|4.2ip is 

(-dzv{k*Vw*'^) =div{Xlu*\K) in K 
\-k*Vw*^'' -nsK ={Ilu*\K)-n9K on dK. 

By a straightforward calculation (ref. [2S), we have the following lemma. 
Lemma 4.8. Letw^ andw*'^ be defined in ^4-21^ and ^4.22\ ), respectively. Then 

—kfVw^ =IhU*\K in K 

(4.23) 

-k*Vw*'^ =IIu*\k in K. 

Moreover, we have 

Lemma 4.9. Letw^^ andw*'^"^ be defined in ^4-21^ and ^.2^ , respectively. Then 

S/wf e xlf^{K), Vw*'"" e Vh{K). 

Proof. It is obvious that 

-div{k,{x)\/wf) = div{Ilu*\K) = ( E /"* ■ nds{divxf)), 

and 



-fce(a;)Vwf • Tie = {^h\K)u* • rie = ( / u* ■ ndsx^-) ■ n^e = { u* ■ nds) 



By Theorem 12.11 and equation p.ip . it follows 



eCdK 

Similarly, we can show that by (|4.20p 



Vu;*^^ = (r)-i {i[u*-nds)x^). 



eCdK 



This completes the proof. □ 

For the relationship between and , we have the following lemma. 
Lemma 4.10. Let be defined in 14.21^ . Then 



\Pe - w'^li.K < Cie + h+ V7h)\\p*\\2,K + CV7h^\\p*\\i^oo^K. (4.24) 

The proof of Lemma 14.101 can be obtained by using Lemma 14.81 and the proof of 
Theorem 3.1 in [14]. 

We have the following convergence theorem on the periodic case. 

Theorem 4.11. Let [{0e,Ue),Pe) be the solution of 112. 3\) and {9'f^,u',^),p'-f^, X'f^) be 
the solution of i3.5\) . Then 

< + Cie + h+ Vdi)\\p*\\2S2 + CW-||p*||i,oo,n. 
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Proof. By Theorem 14. 4[ we can use problem p.3p to perform an error analysis 
for ((0^, . Due to Theorem 14.21 it suffices to choose {{£,h,Vh),qh)) S x 

such that the right-hand side of (j4.8|) is small. 

Set qii ~ VhPcj the projection of onto Q^. Then Poincarc-Friedrichs in- 
equality implies 

\\Pe ~ qhlWn < Ch\p,\i^n- (4.26) 
We choose ^hIk = Vwf^, and the ^/j g X{ ^ by Lemma Due to Lemma 14.101 

\\e. - ailo,n = {J2\P^- ^^lUy <Cie + h + v^)|b*||2,n + C^m,^{A..27) 

Let Vh ■= 1-hU* ■ Note that div{IhU*)\K = {!) Ki the integral average of / over K. 
Then 

\\div{u, - VH)\Wn = (E 11/ - (/)aH1o,/^)^ < Ch\J\,^n. (4.28) 

K 

Moreover, by Lemma ITSl and (|4.24p 



(4.29) 



\\Ue - Vh\\a,K = \\k(Vpe - kfVwf\o,K 

<C{e + h+ ^)\\p*h^K + C^/d^\\p*\\l^^,K. 
Consequently, it follows immediately that 

Ik. - «/J|o,o < C{e + h + V7h)\\p*\\2^n + Cy||b*||i,oo.a. (4.30) 
Combining g^T]), and gSHl), it follows 

||((^,,7.,),p,)-((0^,<),pL)llxxQ <C;i||/||i^o+C(e+/i+x/^)||p*||2,j2+Cyi|b*|^^ 

(4.31) 

Next we employ the technique used in [8] to estimate \\VdPc — A'J|_i /j. Let 
Th := leKT'aP, - X'f,)^l e Then 

f/i • TT-e = 'PaPe ~ ^/t ^ and otherwise. 

By (|4J6l) . it follows 

||f;,|lo,A' + /ix|Miw^/,||o,K < C/illl^aPe - \i\\o,e. (4.32) 
Define r/i = f/j in K and t/j = in il\ K. By the second equation in p.Sp . we have 

(0^, f^)if + (pj,, diwf^)^ = {VoPe - Xl Xh)e. (4.33) 
Since 0^ = Vp^, Green's formula gives 

{Oe,fh)K + {pe,divfh)K = {VoPe " Xh,Pe)e- (4.34) 
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By using and (032), we get 

= - 0[,Th)K + iPe-Ph,divfh)K 

. , (4.35) 

< ll^e - 6'/il|o,A'||'?/i||o,/f + \\Pe - Ph\\o.K\\divTh\\o,K 

< C{h-^\\p, - pi Wo.K + hi\\9e - 9'f,]\o,K)\\V9Pe ' ^hlke, 

which gives 

\\VaP. - XilWe < C{h--^\\p, - p^llo.A- + h'^\\9, - 9],\\o.k). 

Consequently, 

\\T^0Pe^Xi\tij,= J2 E '^K\\V9Pe~Xi\\le 
KeleCdK 

< C( 5: h~^'hK\\p. pUHk + E ''k\\9. 9X,k) (4-36) 

K K 

<C\\p,-p'h\\ln + Ch^9,-9'^,\\l^. 
Owing to Theorem 4.6 in [5D], we have 

\\Pe~pih,n<C{\\{9,,u,)~{elui)\\x+ inf \\p,-qhh,n]. (4.37) 

Combining (|436l) . (|437| and (|43T1) gives 



rap. - A;J1_i,;, < Ch\\,f\\,,n + C(e + /i + Ve;i)|b*||2,n + C^-||p*||i,oo.n(4.38) 

This completes the proof. □ 

Corollary 4.12. Let {9e,Ue) be the solution of 112. 3\) and (6*^, u^) be the solution 
of /TOP. T/ien 



m,u,) - i9lui)\\x < + C{e + h+Veh)\\p*\\2A2 + C\j -\\p*h,oo^n 



Proof. Due to Theorem 4.5 in [20], it follows that 

m,u,)~i9lu'^)\\x<C{ inf \m,u,)-{^^,v^)\\x}. 

Then using the proof of (|4.25p completes the proof immediately. □ 

4.3.2. Analysis of expanded mixed MsFEM using oversampling tech- 
niques. By Theorem 14.111 the local expanded mixed MsFEM without using over- 
sampling techniques causes a resonance error O(y^). By oversampling technique, 
the resonance error 0(^/^) in Theorem (|4.1ip will reduce to 0(-|). In this subsection 
we sketch the analysis for the oversampling expanded mixed MsFEM. 

Let S* D A' be the local domain described in equation (|3.6p . We assume that 
dist{dK, dS) w Hk for the analysis. Let constant cji be defined in p.8p . The constants 
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Cji allow US to introduce an extension operator to S from K, e.g., : Vh{K) — > 
Vh{S) is defined by 

i 3 I 

Similarly, E'^ : X°-^hi^<) — > Xl,^{S) is defined by 

^^(E/5A")-EE/3.-.'^x. 

i i I 

Since dist{dK,dS) ~ Hk, the extension operator E^ has the property 

\\E^{v'')\\o,s~\\v''\\o,K, for Vz;^eV^„(i^) or X^%iK). 

In order to analyze convergence, we define (x) be the solution (up to a constant) 
of the following equation 

(^dtv{h{x)\/wf) =div{ES{X*u*\K)) in S 

\-h{x)Vwf -nss =E'^{I*u*\K)-n9s on dS. ^' ' 

Here we note that IJ-^u*\k is the local interpolation on Vh{K) (i.e., RTo{K)). We 
define the multiscale interpolation operator I^'' : [H-^{il)]'^ — > ^Th fo"^ velocity by 



Then by using a procedure similar to Lemma 14.81 Lemma 14.91 and straightforward 
calculations, we obtain the following lemma. 

Lemma 4.13. Let wf be defined in J^.gffp . Then 

-k^Vw^ = E^'{I^''u*\k) on S, -fc.Vwf Ik on K. (4.40) 

Moreover, Vwf |if — X^ecSif^/e ^* ' ''^ed,s)fi^ g X°\{K), where fj^ is defined in 

We have the convergence result for the oversampling expanded mixed MsFEM. 
Theorem 4.14. Let {{e„u^),p,) be the solution of ^KW and {01' ,ul\pl' 
be the solution of iS.ll]) . Then 

E \\{iOe,U,),Pe) - ((C«r),Pr)lU(K)xQ(K) + WVap, - Kl-l,h 

K (4.41) 

< C{h + e)(!l/lli,o + Wvhs^ + + x/i)(|b*||i,oo,o + ll/llo.o). 

The proof of Theorem 14.141 is presented in Appendix [A] 

4.4. Error analysis in G convergence. In subsection l4.3[ we have investigated 
the case when fcg in (|2.ip is e periodic. However, the expanded mixed MsFEMs 
can also be applied to non-periodic separable scales. In this section we will discuss 
the convergence for the case of separable scales described in G— convergence (ref. 
[27j). G— convergence is more general than the periodic homogenization described in 
Subsection [ 
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A sequence of matrices is G— convergent to k* if for any open set uj d Q and 
any right-hand side / € H~^{uj) in (|2.ip . if the sequence of the solutions in (|2.ip 
satisfies 



■ p* weakly in H^{oj) as e ^ 0, 



where p* is the solution of the equation (j4.15p . in which k* is the homogenized matrix 
in the sense of G— convergence. The G— convergence implies that 

fccVpe fc*Vp* weakly in -/j^(w) as e — > 0. 

There is no explicit formula for the matrix k* , which is defined as a limit in the 
distributional sense, i.e., 

fceVAC fc*e, in 2?'(w;R''), 
where the auxiliary functions Ml (i = 1, • • • , d) satisfy 

Ml ^ Xt weakly in H^{uj) as e 0. 
The auxiliary functions are not explicit. We define the corrector matrix VAC = 

The following lemma is about corrector in G— convergence. 

Lemma 4.15. 1291 Let he a sequence G— converging to k* as e — ^ 0. Then 

Vpe = VAC-Vp*+i?^, 

where — > strongly in L^{oj) as e — > 0. Moreover, ifVM is bounded in L^{io) for 
some r such that 2 < r < oo, and Vp* G L'^i'^) for some s such that 2 < s < oo, then 
i?" — > strongly in L*{uj), as e 0, where t — min{2, ^q^}- We have the following 
convergence theorem for the case of G— convergence. 

Theorem 4.16. Let ((0e,Ue),Pc) be the solution of \2.3\) and {9'-f^,u'-f^,p';^, X'',^) be 
the solution of (3^. IfVNe € L°°{K) for all K, then 

Ihn - {{elui),pi)\\x>cQ + WVoPe - ALlUi,,.} 

<Ch{\u*\o.n + \\f\\i.n). 

Proof. The proof of (|4.42p is similar to the proof of (|4.25l) . We set = VhPe, 
£,h\K = Vuif and Vh ~ IhU* and then use Theorem 14.21 By the proof of (|4.25p . it 
immediately follows 

\\div{u, - Vh)\\o.n + \\Pc- qhhs-i < C/i(b,|i,o + \f\i.n) < Ch\\f\\i,n- (4.43) 

By Lemma [4.81 and Lemma [4. 151 

||Vp, - Vwf ||o,if 

< II Vp, - WK^p*\\o.K + llVAA,(Vp* - Vu;*-^)||o,K + ||VACVu;*^^ - Viff |k.| 

< II Vp, - VAA,Vp*||o,K + II VAC -!,>*)) llo^A' + ||i?f||o,K 

< ||Vp, - VACVp*||o,K + Ch\u*\o.K + WRfWo.K, 
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which gives 

l|Vp, - Vailo^n < C\\R'^\\o,n + Ch\u*\o,n + ||o,k- (4.45) 

K 

Similarly we have 

II". - Vh\\o,n < C\\R%^n + Ch\u*\„,n + C^] \\Rf\kK- (4.46) 

K 

By combining (|4.43p . (|4.45p . (|4.46p and Lemma li?T51 and letting e -> 0, we have 

lim||((0„u,),Pe)-((0L«L),pUllxxQ<C/i(|u*|o,n + ||/||i,o). (4.47) 

By following the proof of Theorem 14.111 and using (j4.47p . we can immediately obtain 
that 

hm WVop, - X[\\_, ,^ < Ch{\u*\o.n + ||/||i,n). (4.48) 

The proof of (|4.42p is completed by combining (|4.47p and (|4.48p . □ 

Because no explicit formula is available for homogenization matrix k* , we are not 
able to obtain an explicit convergence rate in G— convergence. The approximation 
rate is only presented in terms of limit as e — > 0. We can also obtain the conver- 
gence of oversampling expanded mixed MsFEM in terms of G— convergence, and the 
convergence result is the same as in Theorem 14. 161 

4.5. Error analysis for global expanded mixed MsFEM. In this subsec- 
tion, we present the convergence results for the expanded mixed MsFEM using global 
information. We consider the continuum scales for the global expanded mixed Ms- 
FEM. 

For analysis, we make the following assumption for global information. 
Assumption 4.1. There exist global fields {ui,- • ■ ,un} and continuous functions 
Ai{x) (i ~ 1, ■ ■ ■ N ) such that 

N 

u = J2Mx)ui, Ai{x)eC". 

i=l 

The existence of these global fields has been discussed in the literature recently (see 
e.g., [2 [TUl HH ISO]). In particular, Owhadi and Zhang [30] develop global fields 
(iti , . . . , Ud) that are solutions of the following equations 



(4.49) 



where x ~ {xi, . . . ,Xd)- Using these global fields, and Theorem 1.3 in |30j . we can 
prove the following proposition. 

Proposition 4.17. Let {9,u) he the solution of i2.S\) and [Oi,Ui,pi) solve equa- 
tion 14-4S^ for i = 1, . . . , d. Let P = (pi, . . . ,Pd)- //(VP)"^A;VP satisfies the anisotropic 
condition described in Theorem 1.3 in ESI/, then 




^A,{x)u,, e^^A,{x)e^, 



i=i 
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where Ai(x) are continuous functions. Moreover, for any s > d, i = 1, . . . , 

where the constant C depends on the anisotropic condition and s. 

Remark 4.1. To simplify computation of global fields and the coupled coarse- 
scale system by expanded mixed MsFEMs, we may use a single global field (N = 1) 
to construct the basis functions. This simplification has been shown to be effective for 
relevant two-phase flow simulations |IJ 126}/ . 

Let T^u be the interpolation of u, with its restriction on K defined by 

N 



i=l eCdK 

where the velocity basis function ipf^^ is defined in p.l2p . Then it follows that 

{div{u - T^u), Qh) = for Vq/i G Qh- 
By Theorem 3.4 in [1], there exists < /.t < 1, which depends on a and d, such that 

\\u-Ilu\\os,<Ch^{Y,\\A;\\c^,^n))- (4.50) 

i 

Let ^/{x) be the solution (up to a constant) of the following equation 

{—div{kyw) ~ div{If^u\K) in K 
-k\/w ■ HQK = i^h^lx) ■ nQK on dK. 

Then we have the following lemma. 

Lemma 4.18. Let w be solution of Then 

—k{x)\'w=T^u\K on K. 

Moreover, Vw G Xf f^{K). 

The proof of Lemma 14.181 is given in Appendix |Bj By using Lemma 14.181 the 
interpolation estimate (|4.50p and the techniques in the proof of Theorem 14. Ill we can 
derive the convergence result for the global expanded mixed MsFEM. 

Theorem 4.19. Let {{0^,u^),pe) be the solution of ifO|] and {Ol,ul,pl,Xl) be 
the solution of 113.14^ . Under Assumption \4-. 1[ there exists ^ (t) < < 1 j such that 

\\{{e,u),p) - {{elul),pl)\\x.Q + \\Vap-\i\\_.^,, < Ch^, (4.52) 

whereC ^CiWfh^n). 

Theorem 14.191 shows that the convergence of global expanded mixed MsFEM 
is independent of small scales, and the resonance error is removed using the global 
information. 

5. Numerical results. In this section, we present numerical results that high- 
light the advantages of the proposed expanded mixed MsFEMs. In particular, we 
show the improvement in the accuracy of the gradient unknown 9 and the velocity 
unknown u obtained with this new family of methods. In addition, we confirm that 
if the media has strong non-local features (e.g., channels and fractures), then the use 
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of global information in the construction of the basis functions significantly improves 
the accuracy of the multiscale solution 

In all of the numerical experiments reported below, we use RTq in three dimen- 
sions, both as the fine-scale velocity space and as fine-scale shape functions for the 
space approximating Vp in L^{^1)^ (not conforming in H{div,ft)) . The domain ft 
is discretized by a uniform hexahedral fine mesh of size h and a uniform hexahcdral 
coarse mesh of size H > h. The fine mesh is nested in the coarse mesh. The refer- 
ence solutions for pressure, velocity, and Vp are computed by solving the fine-scale 
expanded MFEM. 

Our primary focus in this section is on the behavior of the error of the MsFEM 
solutions for the velocity and pressure gradient. These errors are significantly affected 
by features of the multiscale models and by the different boundary conditions that are 
used in the construction of the multiscale basis functions. The velocity is particularly 
important and required for multi-phase simulations. All errors reported in this section 
are measured in L^. We are also interested in the practical application of simulating 
multi-phase flow. Therefore, in addition to velocity error data, we present solutions 
obtained by the IMPES (implicit pressure, explicit saturation) method for the coupled 
saturation and pressure equations of two-phase flow (incompressible and immiscible). 

5.1. Low permeability channel. In the flrst experiment, the domain ft is the 
unit cube. We take h ~ 1/24 and H — 1/8, so that each coarse element contains 
3x3x3 fine elements. The source function is as follows: 



This source function represents injection and production wells in opposite corners of 
the domain, occupying coarse elements. It is a practical example that demonstrates 
corner to corner flow across the domain, allowing us to test different numerical meth- 
ods in simulating flow for various permeability fields. 

In the case that the permeability is positive and bounded well above zero every- 
where, the expanded mixed MsFEM gives comparable (almost identical) results to 
the standard mixed MsFEM. However, if the permeability is not positive everywhere, 
the standard mixed MsFEM is not applicable. Moreover, if the permeability is not 
bounded well above zero everywhere, the finite precision may significantly degrade 
the its performance. Table \57U demonstrates this fact, where a test is reported using a 
channel of low permeability. This model problem does not have well separated scales. 
The channel has one fine cell's width in the y— and z— directions, and it extends 
across the entire domain in the a;— direction. More specifically, the channel occupies 
(0, 1) X (^, jg) X (^, j^). Inside the channel, k ~ 10~^, and elsewhere k — 1. 

The label "MsFEM" in Table O represents the standard mixed MsFEM, with 
variables for pressure and velocity only. This method makes no direct approximation 
of the pressure gradient, but we approximate Vp with velocity divided by perme- 
ability, when the permeability is positive everywhere. "Expanded MsFEM" is the 
expanded mixed MsFEM, which approximates the pressure gradient Vp locally by 
non-conforming Raviart-Thomas elements. "GMsFEM" means that global informa- 
tion, specifically the fine-scale velocity obtained from the expanded mixed FEM, was 
used in constructing the coarse multiscale basis functions. We observe that using 
global information reduces the error by more than a factor of 20 in both the veloc- 
ity and the pressure gradient, while the difference between standard MsFEM and 
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Table 5.1 
Low permeability channel 





MsFEM 


GMsFEM 


Expanded MsFEM 


Expanded GMsFEM 


Velocity error 


1.173 X 10"^ 


5.005 X 10-^ 


1.173 X 10"^ 


5.005 X 10--'' 


Velocity norm 


5.649 X 10^^ 


5.649 X 10"^ 


5.649 X 10"-^ 


5.649 X 10"'* 


Vp error 


1.739 X 10-^ 


5.284 X 10-^ 


1.858 X 10-'^ 


5.7812 X 10-^ 


Vp norm 


5.670 X 10-^ 


5.670 X 10"^ 


5.670 X 10"^ 


5.670 X 10-^ 



Table 5.2 
Channel with vanishing permeability 





Expanded MsFEM 


Expanded GMsFEM 


Velocity error 


1.140 X 10"^ 


1.122 X 10-^ 


Velocity norm 


5.642 X 10"^ 


5.642 X 10"^ 


Vp error 


1.461 X 10"^ 


9.338 X 10-* 


Vp norm 


5.782 X 10-^ 


5.782 X 10-^ 



expanded MsFEM is negligible. In the remainder of this section, we do not report 
results for standard MsFEM. 

Next, we test the more interesting case of a permeability that vanishes in the 
channel. Standard mixed MsFEMs cannot be used in this case. In each cell within 
the channel, we define the permeability to be the vanishing, non-negative function 

A: = max(0, l-32a;(l-.T)2/(l-2/)), for (x, y) € (0, 1)^ (5.2) 

defined on the reference unit square (see Figure [01 for the contour plot). Outside the 
channel, fc = 1. This permeability has no scale separation, and it ensures a positive 
definite mass matrix (weighted by fc) for the discontinuous space approximating Vp 
in L^{rt)^. Note that expanded mixed MsFEMs allow the permeability to vanish 
locally within fine cells, unlike standard mixed MsFEMs which involve k~^, but the 
mass matrix weighted by k must be invertible. The results for the expanded methods 
are shown in Table 15.21 Global information reduces the relative velocity error from 
approximately 20% to 2% and the relative pressure gradient error from approximately 
25% to 16%. 

5.2. Oscillatory permeability. The next experiment demonstrates the reso- 
nance error that appears in the velocity and the pressure gradient that are obtained 
using the local expanded MsFEM. The permeability is taken to be the oscillatory 
function 

k = (sin(207rx) + 1.5)(sin(207r2/) + 1.5) (5.3) 

on the domain fl = (0, 1) x (0, 1) x (0, i), and the mesh has 100 x 100 x 8 fine cells. The 
source function is the two-dimensional analog of (j5.ip (constant in the z-direction), 
namely 



r 1 if (a;,y) isin (0,i7)2, 
fix,y,z)=< -1 if (a;,y) is in (1 -iJ, 1)2, 
[ otherwise. 



(5.4) 
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Fig. 5.1. Partially vanishing permeability within a reference fine cell 

Table 5.3 

Oscillatory permeability without scale separation, 100 X 100 x 8 fine mesh 





20 X 20 X 1 coarse mesh 


10 X 10 X 1 coarse mesh 




Exp. MsFEM 


Exp. GMsFEM 


Exp. MsFEM 


Exp. GMsFEM 


Velocity error 


2.066 X 10-^ 


1.050 X lO-'' 


1.827 X 10-^ 


7.606 X 10-* 


Velocity norm 


5.143 X 10'^ 


5.143 X 10"^ 


5.143 X 10"^ 


5.143 X 10-^ 


Vp error 


1.370 X 10-^ 


1.452 X 10"" 


1.552 X 10"^ 


5.360 X 10-* 


Vp norm 


3.254 X 10-^ 


3.254 X 10-^ 


3.254 X 10-^ 


3.254 X 10-^ 



Thus the problem is essentially two-dimensional, although we solve in a three-dimensional 
domain. The solution will be approximately constant in the z-direction, so we simply 
plot an x-y slice. We tested both local and global expanded mixed MsFEMs with 
coarse meshes of size 20 x 20 x 1 (each coarse cell having 5x5x8 fine cells) and 
10 X 10 X 1 (each coarse cell having 10 x 10 x 8 fine cells). In the 20 x 20 x 1 case, 
the permeability has one period per coarse edge length. Table 15.31 lists the error for 
both methods on the two choices of coarse grids. 

The local expanded mixed MsFEM clearly suffers from the resonance errors in 
the velocity and the pressure gradient. In this example, refinement of the coarse mesh 
does not change the error significantly because the oscillations in k are not resolved, 
even by the 20 x 20 x 1 coarse mesh. In fact, the error of the velocity actually 
increases slightly when the coarse grid is refined. Figure 15.21 shows the a;-component 
of the reference velocity on the fine grid (top left), and the excellent agreement with 
the x-component of the velocity obtained with the global expanded mixed MsFEM on 
the 20 X 20 X 1 coarse grid (top right). In contrast, the impact of the resonance errors 
on the x-component of the velocity obtained with the local expanded mixed MsFEM 
is evident in Figure [5?2] for both the 10 x 10 x 1 (bottom-left) and 20 x 20 x 1 (bottom- 
right) coarse grids. Note that these plots are in an x-y slice, as the solutions are 
constant in the z-direction. On the other hand, the global expanded mixed MsFEM 
performs well in this example, and the refinement of the coarse mesh reduces the error 
by more than half. 

5.3. IMPES with random shales. To compare the performance of the local 
and global expanded mixed MsFEMs for a two-phase fiow problem, we use an Im- 
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Fig. 5.2. The x-component of the velocity for the oscillatory permeability field 1 15.311 is plotted 
for the fine- scale reference solution (top-left), and for three expanded mixed MsFEMs. The solution 
obtained with the global expanded mixed MsFEM on the 20 X 20 X 1 coarse grid (top right) shows 
excellent agreement with the reference solution. The solutions obtained with the local expanded mixed 
MsFEM on the 10 X 10 X 1 coarse grid (bottom left) and the 20 X 20 X 1 coarse grid (bottom right) 
show the impact of the resonance errors. 



plicit Pressure Explicit Saturation (IMPES) formulation. The boundary conditions 
that we use to define the global expanded MsFEM basis functions are obtained from 
the fine-scale expanded MFEM solution at the initial time. We note that in very 
complex highly varying flows it may be advantageous to update these basis functions 
by updating the fine-scale expanded MFEM solution that is used in their construc- 
tion. However, our current study is focused on the effect of the two different velocity 
approximations, and hence, we fix the multiscale bases throughout the simulation. 

The two-phase fiow scenario that we consider is the traditional two-spot problem 
on a square domain (e.g., [U [25]), in which the water is injected at the bottom left 
corner and oil is produced at the top right corner. Here, we use simple fixed rate 
wells that are constant over a coarse cell. We note that it is possible to define a well 
on a fine cell, and capture the infiuence of the near-well heterogeneity, by modifying 
the source term in the construction of the basis functions. Similar modifications were 
developed to address this issue in the subgrid upscaling method [4] and the standard 
mixed MsFEM [2]. However, fine-scale wells arc beyond the scope of this paper. 

This two-spot problem is essentially two dimensional, with the fine-grid discretiza- 
tion having h = 1/100 (100 x 100 x 8 mesh) and H = 1/10 or 1/20. The logarithm of 
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Fig. 5.3. Logarithm of permeability field used in the two-phase flow two-spot problem. 
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Fig. 5.4. Relative saturation errors (percentage) at each timestep of the IMPES simulation, 
for global and local expanded mixed MsFEM 

the permeability field is plotted in Figure [Ol and includes features that are represen- 
tative of shale barriers and sandy deposits or channels. Specifically, the darkest blue 
cells have a permeability of 10~^, but are masked by the vanishing function given in 
(|5.2p on the reference cell to represent a shale barrier. In contrast, the red cells have 
a permeability of 100 and represent sandy deposits or channels. The remainder of 
the finc-scalc cells are assigned a moderate or low permeability, and are shown as two 
distinct colors. 

The water and oil mobilities are defined simply as the quadratic functions Xw{Sw) = 

— and Xo{Sw) = respectively. Here, denotes the viscosity of phase a, 

and we take fiw = 0.5, fio = 1- The total mobility A(S'„,) := A,u,(S'„,) + \o{Sw)- Then 
the pressure equation is given (in the absence of gravity and capillary effects) by 

—div{\{Stjo)k\/p) = q. 

At each timestep, we use the velocity u computed from the pressure equation to ad- 
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Fig. 5.5. Water saturation at the 700th timestep computed with the fine-scale reference velocity 
(left); with the coarse global expanded mixed MsFEM velocity (middle); and with the coarse local 
expanded mixed MsFEM velocity (right). 



Table 5.4 

Randomly vanishing permeability without scale separation, 100 X 100 X 8 fine mesh 





20 X 20 X 1 coarse mesh 


10 X 10 X 1 coarse mesh 




Exp. MsFEM 


Exp. GMsFEM 


Exp. MsFEM 


Exp. GMsFEM 


Velocity error 


2.184 X 10-^ 


5.438 X 10-'' 


2.218 X 10-^ 


7.137 X 10-* 


Velocity norm 


6.202 X 10-^ 


6.202 X 10"^ 


6.202 X 10-^ 


6.202 X 10-^ 


Vp error 


3.529 X 10"^ 


1.999 X 10"^ 


3.636 X 10"^ 


2.283 X 10-^ 


Vp norm 


1.471 X 10"^ 


1.471 X 10-^ 


1.471 X 10-^ 


1.471 X 10-^ 



vance the water saturation Sw by an explicit discretization of the saturation equation 
^^"^ + div{f^{Sn,)u) = 0, with fwiSw) ~ ^wiSw) 



Figure [5^ shows a plot of the saturation error in the L^(il) norm (relative to the 
reference IMPES saturation computed with the fine-scale expanded mixed FEM veloc- 
ity), at all 2000 time steps. The water saturation, Sw, at the 700th timestep is plotted 
in Figure [575] bv using the fine-scale velocity and the global and local expanded mixed 
MsFEM velocities. Figure 15.51 shows that the IMPES solution based on the global 
expanded mixed MsFEM velocity remains close to the reference solution throughout 
the long simulation, whereas the local expanded mixed MsFEM velocity results in 
substantial error very early in the simulation. The error naturally decreases later in 
the simulation, when most of the domain becomes highly saturated with water. 

The IMPES experiments discussed above used coarse velocity fields on a 10 x 10 x 1 
coarse mesh (H = 1/10), with each coarse cell discretized by 10 x 10 x 8 fine cells. 
To study the coarse velocity error behavior, we also computed the coarse velocity on 
a 20 X 20 X 1 coarse mesh, with each coarse cell discretized by 5 x 5 x 8 fine cells. 
The results are reported in Table 15.41 The global expanded mixed MsFEM velocity 
is about 25% more accurate on the 20 x 20 x 1 coarse mesh than on the 10 x 10 x 1, 
while the local expanded mixed MsFEM velocity improves by less than 2%. Without 
global information, non-local effects result in large errors even after refinement. 

We conclude that the expanded mixed MsFEM method is very effective in ap- 
proximating both the pressure gradient and velocity when a locally vanishing perme- 
ability field makes standard mixed MsFEM methods infeasible. Practical simulations 
of two-phase flow demonstrate the high degree of accuracy attained by using global 
information in constructing the coarse basis functions. Thus the expanded mixed for- 
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mulation with global information is the most robust and accurate method considered 
in this study for applications involving a locally vanishing permeability field. 

6. Conclusions. We developed a family of expanded mixed MsFEMs for elliptic 
equations and considered their hybrid formulation. In this formulation, the four 
unknowns were solved simultaneously, namely, the pressure, gradient of pressure, 
Lagrange multipliers, and velocity. The expanded mixed MsFEMs work well for 
the case that the coefficient (e.g., permeability) is very small or locally vanishing 
in some regions of the domain. This case has important applications in reservoir 
simulation. In contrast the standard mixed MsFEMs require the coefficient to be 
uniformly positive everywhere. In addition, expanded mixed MsFEMs provide the 
gradient of pressure without significantly increasing the computational cost relative 
to the standard method. 

We analyzed the expanded mixed MsFEMs for separable scales and non-separable 
scales, and established a priori error estimates. We showed that the convergence 
rates of the local expanded MsFEMs depend on both the small physical scales and 
on the coarse-mesh size. In contrast, we showed that the global expanded mixed 
MsFEMs achieves a convergence rate that only depends on the coarse-mesh size. To 
support this analysis we applied the expanded mixed MsFEMs to various models of 
flow in heterogeneous porous media. For the models with non-separable spatial scales, 
significantly better accuracy was achieved by using global information to construct the 
multiscale basis functions. Moreover, for the periodic example of separable scales using 
global information eliminated the resonance errors observed in the the local expanded 
mixed MsFEM. These numerical results showed the efficiency and robustness of the 
developed methods. 

Acknowledgments. We would like to thank the referees for their comments and 
suggestions to improve the paper. 

Appendix A. proof of Theorem I4.14i 

Due to Theorem 14.51 problem p.ll|) is equivalent to problem p.lOp , and we can 
use problem p.lOp to estimate the errors for {{Of^ TU°f^),p1^). By the estimate (|4.9|) . 
we have 



<C'{,, inf ^ ^^\\{{0,,u,),p^)-{{£,h,Vh),qh)\\x{K):<Q{K)] 

^ (A.l) 



E 



= h+l2- 



J2k \\''^h\\div,K 



We set Qh = VhPe, ^h\K = Vwf lif and Vh = ^h^u* in (|A.1[) . Lemma 14.131 shows that 
£,h G ^1 ft ^^'^ '^h G ^Th- ^^^^ Lemma 4.6 in [14], we obtain 

|/i| < C{h + e)(||/||i,f2 + \\p*h.n) + C(J + x/i)||p*||i,oo,a. (A.2) 

By using the proof of Theorem 2.2 in [T3], we estimate the consistence error I2 by 

I/2I < C{h + e){\\f\\,s, + \\p*hs{) + + V^)(!|p*||i,oo,o + ||/||o,o). (A.3) 
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Next we estimate \\V9p, - A^,*||_i_^. Let th := \e\{VoPe - K')^x. ^ ^2%W and 
the operator A^/i be defined in p.9p . Then 

■MhiTh) ■ n-e ~ VdPf: — A^* on e and otherwise. 

By direct calculation or scaling argument, it follows 

\\Mh{rh)h.K + hK\\divMh{Th)\WK < ChlWVsp, - Ari|o,e. (A.4) 
Define Th = Th in K and t/j = in 0\ K. By the second equation in p. lip , we have 

{ei\ Th)K + {pl'AivTh)K = {VaP. - Ar , Ar )e. 
Since divM.h(jh) = divfh, it follows that 

{9r,rh)K + {pl',divMh{Th))K = {Tap. - K, Ar )e. (A.5) 
Since 0^ = Vpc, Green's function gives 

{e,,Mh{Th))K + {pe,divMh{Th))K = {VoPe ~ K^Pe)e. (A.6) 

By a similar argument to (|4.6p . it follows that 

\\vh\WK<C\\MhVhh.K for -ivheX^^iK). (A.7) 
By using (|A.5p . (|A.6I) . (jA.4p and (|A.7p . we get 

rap. - Arllt = iT'dP^ - K'^-POP. - \V)e = (Pe - Xl^VoPe " A^e 

= {6, - 9l\Th)K + {9„Mh{Th) - Th)K + {p,-pl%divMh{Th))K 

<\\0.- ei'\\o,K\\^hh.K + \\e,\WK{\\Mh{Th)\WK + \\Thh,K) 

+ \\p,-pl'\\o,K\\d^vMh{fh)\WK ^^'^^ 

< C\\e, - 0rilojdl>'/.(^/OI|o,X + C\\e,\\o^K\\Mh{fh)h.K + \\pe~pl'\\o,K\\dlvMh{fh)\\o,K 

< C{h~J\\p,^pl'\\o,K + hi\\0, - CI|0,A- + hi\\e,\\o,K)\\VaPe - Arilo.e, 

which gives 

WVaPe - Ar ||0,e < C{h-^\\p._ - pnO,K + 4ll^. - dl0,A' + 4ll^e||0,i^). 

Consequently, 



ll^oPe-Arll^,, = E E hK\\Vap.-Xl^ 

Ke%h eecdK (A.9) 



\ OS 1 1 2 
l0,e 



< c\\p, - pnln + ch^e, - + ch^e.WoM- 

Combining (|XT1) . (|X2l) . (|X3l) and ((X9)) proves dMl]). 

Appendix B. proof of Theorem 14.181 

By straightforward calculation, it follows that 



fc'«"Ej=iEecOK(/e^'"' • '^'=ds)(-fcV(/)f^^J] =div{I^u\K) mK 
lE^iEecaxl/e^'"' •"erfs)(-fcV<?!)f^J -n =(I^it|K)-n on aA'. 
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Comparing equation (|4.5ip with equation (|B.1|) . we find that 

N 

-k\/w = ^^{ A,u, ■ neds){-k\/(t)^^J = IIu\k. 

i=l eCdK •'^ 

By the uniqueness of solution for (|4.5ip . it follows that 

N 

i=l eCdK 
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